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ABSTRACT 

The theory and application of numerical methods for unstructured meshes 
have been improved significantly in recent years. Because the grids can be 
place arbitrarily in space, unstructured meshes can provide much higher spatial 
resolution than regular meshes. The built-in nature of mesh adaptivity for 
unstructured meshes gives one way to simulate highly dynamic, hierarchical 
problems involving both coUisionless dark matter and collisional gas dynamics. 
In this paper, we describe algorithms to construct unstructured meshes from a 
set of points with periodic boundary conditions through Delaunay triangulation, 
and algorithms to solve hydrodynamic and N-body problems on an unstructured 
mesh. A combination of a local transformation algorithm and the traditional 
Bowyer- Watson algorithm gives an efficient approach to perform Delaunay 
triangulation. A novel algorithm to solve N-body equations of motion on an 
unstructured mesh is described. Poisson's equation is solved using the conjugate 
gradient method. A gas-kinetic scheme based on the BGK model to solve 
Euler equations is used to evolve the hydrodynamic equations. We apply these 
algorithms to solve cosmological settings, which involve both dark and baryonic 
matter. Various cooling and heating processes for primordial baryonic matter 
are included in the code. The numerical results show that the N-body and 
hydrodynamic algorithms based on unstructured meshes with mesh refinement 
are well-suited for hierarchical structure formation problems. 

subject headings numerical methods, cosmology, galaxy formation. 

1. Introduction 

Numerical simulations in astrophysics turn out to be very challenging because of the 
large dynamical range required in three dimensions. Examples include modeling of star 
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forming regions and the origin of galaxies. In cosmology, structures are believed to have 
formed hierarchically, requiring a simultaneous modeling of structures on scales of ~ 100 
Mpc and ~ 10 kpc. Various hydrodynamical techniques have been explored to achieve 
such a large dynamic range, from Eulerian methods using regular meshes (c.f. |Cen 199^ 



and [Ryu et al. 19931) and recently with Adaptive Mesh Refinement (AMR, c.f. [Berger fc 



[Colella 1989| , [Klein, Colella, fc McKee 1992|) to Lagrangian methods like Smoothed Particle 
Hydrodynamics (SPH, c.f. [Hernquist fc Katz 1989| , [Katz, Weinberg, fc Hernquist 1996[) . 



Eulerian schemes without AMR are inadequate because of computational expense and 
are wasteful because high resolution is typically not required at all points in a simulation 
volume. SPH methods provide higher spatial resolution than regular mesh Eulerian methods 
such as TVD and PPM, but have poorer shock resolution than shock capturing methods. 

N-body codes can be classified as direct, in the case of Particle-Particle methods 
and TREE methods arnes fc Hut 198^ , [Hernquist 1987[) , or grid-based, such as the 



Particle- Mesh method (pUfstathiou et al. 1985[) , or hybrids, such as P^M ( [Hockney fc 



Eastwood 1981[ ) or TPM { \Ku 1995[ ), depending on the potential solver. In cosmological 



simulations involving gas and dark matter it is desirable that the N-body and hydro solvers 
achieve similar spatial resolution. Normally, Eulerian schemes for N-body and gas are 
combined, such as PM+TVD code ( [Ryu et al. 1993[ ), and Lagrangian schemes go together, 
such as TREESPH ( Hernquist fc Katz 1989| ). Numerically it is appreciate to employ similar 



algorithms for N-body and gas dynamics. 

Recently, unstructured meshes have become increasingly popular in many fields, such 
as geophysics, seismology, structural mechanics and computational fluid dynamics. When 
combined with an accurate shock-capturing technique, codes employing unstructured 
meshes have many advantages over particle based algorithms and Eulerian codes with 
regular grids, in principle. In an unstructured mesh, grid points are connected by triangles 
in two dimensions and tetrahedra in three dimensions. Since grid points can be placed 
arbitrarily, an optimal mesh can be configured for any applications. Mesh refinement can 
be achieved by simply adding more grid points and reconnecting the mesh. The refined 
mesh will have the same topology as the original one except that it will have more cells, 
thus mesh refinement adds no overhead to algorithms designed for unstructured meshes. 
Although the nodes in an unstructured mesh can be irregularly distributed, the internal 
data structure used to represent the grid is homogeneous as opposed to block-structured 
grids where block boundaries and interiors must be distinguished. As described below, only 
local grid operations are needed to solve equations on unstructured meshes, hence the codes 
can be easily parallelized. 



In this paper, we describe techniques to construct unstructured meshes and algorithms 



- 3 - 



to solve N-body and gas dynamic problems on these grids. We apply this code to 
perform cosmological N-body + hydrodynamic simulations, including various cooling and 
heating processes and mesh refinement. To avoid lengthy mathematical derivations of the 
algorithms, we put all the required formulas in the appendices. 



2. Numerical Algorithms on Unstructured Meshes 

2.1. Construction of Unstructured Meshes 

Unstructured meshes are constructed from arbitrarily scattered points in n-dimensional 
space by Delaunay triangulation. The mathematical definition of Delaunay triangulation 
can be found in text books on geometric design. An example of Delaunay triangulation 
in 2-dimensional space is illustrated in Figure |I|. The scattered points are connected by 
non-overlapping triangles obeying certain rules. The interior of the circumcircle of any 
triangle contains no other point in the point set. Delaunay triangulation is unique provided 
that no + 2 points are co-spherical in n-dimensional space. There are many properties 
associated with Delaunay triangulation (c.f. [Lawson 1986| , |Barth 1995| ), many of which will 



be cited in our paper without strict mathematical description. In Appendix A, we gave the 
essential formulas for geometric relations between a point and a simplex (triangle in 2D and 
tetrahedron in 3D). 

Much effort has been devoted to designing algorithms to perform Delaunay 
triangulation. Among these, incremental insertion algorithms are of particular interest, 
since we will incorporate mesh refinement to enhance resolution. The methods of Bowyer 
Bowyer 1981| and Watson [Watson 1981| are very similar. When a new point is inserted 



into the existing triangulation, those simplices with their circumspheres enclosing the new 
point are deleted and and new simplices corresponding to those just deleted are added. 
The Bowyer- Watson algorithm is straight-forward to implement and is efficient {0{N^~^^^'^) 
in n-dimensional space). Another triangulation method is the edge swapping algorithm of 
Green and Sibson Preen fc Sibson 1977| , which was extended to 3D by Joe [Joe 1989] . Edge 



swapping algorithms insert a new point into the simplex that encloses it, and perform a 
sequence of local transformations until no further local transformations can be performed. 
The edge swapping algorithm is slightly slower (also 0{N^~^^^^)), but can be adopted to 
handle different criteria to perform local transformations. 

For cosmological simulations, we require a triangulation algorithm for a set of points 
in a periodic box. The Bowyer- Watson algorithm appears to have difficulties with periodic 
boundaries when the number of points in the box is small, because some triangles might 
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have two vertices that are periodic images of the same point. The edge swapping algorithm, 
however, can be apphed to periodic volumes. Thus a combined Bowyer- Watson algorithm 
with edge swapping algorithm is efficient for periodic boxes. In Figure |^ we show an 
example of a 2-dimensional unstructured mesh for 1000 random points in a periodic box. 

Our data structures to describe the geometric connections in an unstructured mesh 
differ from those used previously (e.g. [Kallinderis fc Vijayan 1993| ). Two data types 



represent an unstructured mesh: nodes and cells. A node contains a vector to describe the 
position of a point and a flag to record the refinement level. A cell contains n + 1 pointers 
to its vertices, another n + 1 pointer to its neighbor simplices, and an integer flag to record 
various information about the cell's status and its relation to its neighbors. With periodic 
boundaries, another two integer flags record the relative position between the simplex and 
its vertices and its neighbor simplices. So the memory requirement with periodic boundaries 
is n + 1 words for each node and 2n + 5 words for each cell. This implies that the total 
memory required to store an unstructured mesh is about 3 + 2x9 = 21 words per node in 
2-D, and 4 + 6x 11 = 70 words per node in 3-D. 



2.2. N-body algorithms on unstructured meshes 

Gravitational accelerations in N-body systems can be calculated either from particle- 
particle methods or particle-mesh techniques. Below we introduce a new particle-mesh 
algorithm for unstructured grids. 

The discretized Poisson's equation on an unstructured mesh can be derived as follows. 
Consider node as illustrated in Figure 0. We choose the control volume for node with 
the boundary connected by the middle points of the edges and middle point of each simplex 
associated with node 0, such as that indicated by the dotted lines in Figure |^. Integrating 
Poisson's equations over this control volume gives, 

/ dVV^(j) = f dVAnGip - pt) = AnG{mo - ptVo), (1) 
Jn Jn 

where Vq and mo are the control volume and mass of node 0, respectively. If we linearly 
interpolate the potential field in each simplex Tj, we have = J2k'>^k{x)(f)k, where 
Wk{x) is the weighting to vertex k. It can be shown that Wk{x) is equal to the barycentric 
coordinates of a point located at x relative to vertex k of simplex Tj (see Appendix A for 
more details). Using Gauss' theorem, the left hand side of equation (|l|) can be written as. 



/ dvv'(j) = ds-v<p = J2s,-v<p = -J2 -^wo ■ E ^Mx)<Pk, (2) 
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where Vi is the volume of simplex Tj and n is the dimension of space. This result is the 
same as that derived using the Garlerkin finite element approximation ( [Barth 1995| ). 

The discrete Poisson's equation resulting from this procedure is a system of linear 
equations, Mij ■ (pj = bi, where Mij is a symmetric sparse matrix. Since the matrix M is 
sparse, the conjugate gradient method can be used to solve this linear system efficiently. To 
guarantee convergence, we require that M be positive definite. A necessary and sufficient 
condition for M to be positive definite is that the triangulation is a Delaunay triangulation 
( Parth 1995|) . The convergence rate of this simple conjugate gradient method is shown in 
Figure ^. The error decreases exponentially with the number of iterations, given a good 
initial guess of the potential. In N-body simulations, we can always use the solution in the 
previous step as the initial guess. Our numerical experiments indicate that the number 
of iterations required to achieve convergence during the next time step is typically about 
20 - 50. 

In order to solve Poisson's equation on an unstructured mesh for N-body problem, we 
must interpolate particles to the mesh. For regular mesh, particle interpolation can be done 
using the Cloud- in-Cell (CIC) interpolation (see, for example [Efstathiou et al. 1985| ). For 
an unstructured mesh, a similar procedure can be used. For each particle, we identify the 
cell containing the particle, and calculate the barycentric coordinate bi of the particle inside 
the cell according to equation ([A^). Mass is assigned to each node of the tetrahedral cell 
with weighing factor 6j. After particle mass interpolation, the density of each node obtained 
by dividing its mass by its dual volume. 

After solving for the gravitational potential, the acceleration on each node is calculated 
from the average of its control volume. For example, the acceleration at node in Figure ^ 
is, 

-> 1 r V 

Fo = - {-V<P)dV = E 77 E '/'aV^fc. (3) 

The acceleration on each particle is calculated in a similar fashion as in the Particle-Mesh 
algorithm, 

n 

Fi = J2wk{xi)Fk. (4) 

k=0 

The above formulation, when applied to a regular mesh, is identical to the Particle-Mesh 
algorithm with CIC interpolation. 

In Figure |^, we show the force between two equal mass particles with different 
separations and orientation obtained using the above algorithm. The force behaves similarly 
to the PM method, i.e. accurate long distance forces, but underestimate short range force. 
The force is not very noisy within one cell, which indicates the force resolution can be 
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improved using techniques similar to P^M. The force resolution is roughly 1.5 — 2 cells, 
which is slightly better than the PM algorithm because each node is connected with more 
cells in an unstructured mesh than in a regular mesh. The deviation from r~^-law at large 
separations is due to periodic boundaries. 

As discussed later, mesh refinement can be performed at low cost, so we can achieve 
high spatial resolution using unstructured meshes for the N-body problem. 



2.3. Hydrodynamics on unstructured meshes 

The Euler equations are solved on unstructured meshes through the Finite Volume 
scheme (c.f. [Vijayan fc Kallinderis 199^) , 



— rfV^+ / F(U)-dS = 0, (5) 
n at Jan 

where U = {p, pv, E}^ represents the fluid state and F{U) is the flux vector. Before we 
write down the discrete form of the equations, we need to decide whether we want to store 
the fluid variables on the nodes or in the cells. The node representation uses less memory 
because the number of cells is typically 5 — 6 times the number of nodes in 3-D. But there 
are indications that the cell representation gives higher resolution than node representation 
(see [Mavriplis 1992| for a discussion). We choose node representation to minimize memory 




usage. Consider node in Figure ^ with its control volume fl illustrated in the flgure. We 
have, 

= TT X! -^o/c ■ Sok, (6) 

where A/J is the set of neighboring vertices connected with node 0, Fq^ is the flux at the 
middle of edge — k, and 5*0^ is the total surface area of the control volume related to edge 
— k. All we need is a method to calculate the flux in the middle of each edge accurately. 

Previous approaches to solve the Euler equations have used Upwind schemes of various 
types to calculate the flux at the edges (c.f. |Barth 19"95| ). Upwind schemes require artificial 



viscosity to achieve better than first order accuracy. Here we introduce a new approach 
based on the gas kinetic theory of fiuid dynamics. 

The hydrodynamic equations (both the Euler equations and the Navier-Stokes 
equations) can be derived from Boltzmann's equation through the Chapman-Enskog 
procedure (see, for example, |Shu 1991| , Chapters 2 & 3). It is quite physical to derive 
numerical schemes for hydrodynamic equations using the gas-kinetic theory. Recently, 
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& Prendergast (1994) and pCu, Martinelli, fc Jameson (1995)| successfully developed such 
numerical schemes based on the BGK model ( [Bhatnagar, Gross, fc Krook 195^ of the 
coUisional Boltzmann's equation, 



9f 9-f 

Ot T 



(7) 



where g{t,x,u) is the equilibrium distribution function and r is the coUisional time scale. 
The BGK model accurately describes a large range of situations, from very high density, 
high temperature flows to very high Mach number (> 10^) flows. The solution to this 
equation can be written as 

1 rt 



fit, X, u) 



T JO 



9 



it', x', u)e-(*~*')/"rft' + e-*/Vo(5 - ut, u) 



(8) 



where x' = x — u{t — t') and /o(x, u) is the initial state. 



The macroscopic quantities U (t, x) and F{t, x) are moments of the distribution function 

f(t,x,u), 



^afdud^ 
uipafdudC, 



where, 



1 

U 



(9) 
(10) 

(11) 



with ^ representing the internal variable with K degrees of freedom which will be discussed 
in Appendix B. 

Consider an edge which connects two nodes. Without loss of generality, we assume 
that the two ends of the edge have coordinates (—1/2,0,0) and (1/2,0,0), and that the 
current time is t = and the time step is AT. the current position in consideration is 
X = 0. Following the treatment in pCu, Martinelli, fc Jameson (1995)1 , can expand the 
distributions /o and g along the x-direction around the edge center as follows. 



fo{x,u) 
g{t,x,u) 



9^ 



g^il + A^x), X <0 
g^il + A^x), X > 

1 + A^^x + Bt, X <0 
1 + A'^^x + Bt, X > 



(12) 
(13) 



Here go is the equilibrium state, which is a Boltzmann distribution for hydrodynamic 
equations. The coefficients of A^, A^, A'^^, A'^^, B can be expanded in velocity space 
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as = Aj^ipp, = A^ipp, A^^ = A^^ip^, A^^ = A^^ipp, B = Bpipp, with 
Aj^, Aj^, A'^^, A'^^, Bp being constants. The reason for the sphtting of the right hand side 
and the left hand side has its physical basis ( pCu, IVlartinelli, fc Jameson 1995| ) and serves as 
the mechanism for shock capturing in the scheme. 

Substituting the solutions to equation (||) after the integration, we get, 

/(t,0,n) = {l-e-'/^)g^J'' + {-T+{t + T)e~'l^)u,A^g^J'' 

+ {t-r + re-*/-)5^7^/^ + e-''^ h{-ut,u). (14) 

This distribution is used to calculate the flux function Fa{t,x) through moment integration 
(equation |TD|) . 

Notice in the above derivation that it does not matter if the gas is 1-D or 3-D. In this 
sense, the BGK gas-kinetic scheme is truly multiple dimensional without involving vector 
splitting which are usually used in TVD or PPM codes to generalize from 1-D space to 
multiple dimensional space. 



3. Cosmological Equations 



The dynamical equations for dark matter and gas in a comoving frame can be written 
as follows (c.f. Peebles 1980| , |Cen 1992|) , 



dpv 



dt 

dvi a ^ 
at a 



dp_ ^ 1 dpvk 
dt a dxk 



1 ^ 

-Vi 

a 

--V0 

a 

(P-Po) 

Qj 



a dxj 



dt a dt 

dE 1 d{E + P)vk 
'dt^ 



-pv ■V(j) + n-A 

a 



a 

-2-E 

a 



a dxk 

The polytropic equation of state is normally adopted for an adiabatic gas 



E = -pi' 



P 



7 



(15) 
(16) 
(17) 
(18) 
(19) 
(20) 

(21) 
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These equations can be simplified by defining, dt' — a ^dt , v'^^ — av^ , E' — a^E , p' — a^p 
and (f)' — a^0. yielding, 

dpv'^ I dpv-v'k ^ dp' ^ _ d^^ 

dt' dxk dxj dxj ' 

dt' dxk dxk 

E' = + (25) 

z 7 — 1 



3.1. Time Integration Scheme 

The equations of motion for the dark matter are integrated using a time-centered, 
second order accurate leapfrog algorithm. The particle positions are one half time step 
ahead of their velocities. 

£."+V2 ^ £.n-l/2 ^ a-l^-n^i, (26) 



- 1 2"V"/ rn , -1 T^n+i/^ /27\ 



where F = — V0. 



In applying mesh refinement, we allow the system time step based on gravity to adjust 
according to, 

I a?6li ■ 
^ max(|Fi|, |F^-|) 

where 5lij is the length of the edge between nodes i and j in the unstructured mesh, and 
Fi is the gravitational acceleration at node i. The constant Cg^av has a meaning similar to 
the CFL condition in hydrodynamics. Our numerical experiments show that Cg^av ~ 0.3 is 
a good choice. 

When the time step changes from 5ti to 5^2, we adjust the particles positions from 
t + Sti/2 to t + St2/2 using the following second order accurate formula, 

Xi{t + St2/2) = Xi{t + 5ti/2) + ^^^l^^Siit) + M^Mlf .(i _ st,/2). (29) 

Z o 



The Courant-Priedrichs-Lewy (CFL) stabihty criterion determines the hydrodynamic 
time step for the system. We use simplified version of the CFL criterion in an unstructured 
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mesh, 

5ihydro < mill ( j^r^ — \ r I , (30) 

13 \msiXk=ij{Vk ■ riij + Cs,k> J 

where 5lij is the length of edge i — j, rtij is the unit vector indicating the edge direction and 

Cs is the sound speed. We argue that the above criterion is sufficient to satisfy the CFL 



stabihty criterion described in Parth (1995 



The gravitational terms in the cosmological hydrodynamic equations can be solved 
consistently in the gas kinetic scheme by including the force term in Boltzmann's equation. 
But it would be rather expensive to do so. Instead, we treat these terms as source terms due 
to an external force. The fluxes due to gravitational acceleration are calculated as follows, 

A^p = 0, (31) 
A^ptT = l(p(-+i)+pW)F'("+i/2)^^/^ (32) 

A^E' = i(pt;'(«+i)+pt;'W)F'("+i/2)5t', (33) 

where F' = — V0'. Since the hydrodynamic quantities are synchronized with the velocity 
field of the dark matter, when the system time step changes, we still need only to change 
the particle positions. For the hydrodynamic time step, we allow for variable CFL constant 
from one time step to another in order to limit the change of system time step. 



3.2. Radiative Cooling 

Various cooling and heating processes relevant to primordial gas have been included 
in the code (see Appendix D for a list of processes). Since the cooling time can be very 
short compare with hydrodynamic time (c.f. Figure |^), we have to be very careful with 
time integration of the energy equation when cooling processes are turned on. In our 
implementation, we integrated the cooling function with adjustable time steps within one 
system time step. The variable step, fifth order accurate Runge-Kutta integrator described 
in Press et al. (1992"]| is used to integrate the following equation, 

du Au ^ , . , 

where u is the thermal energy. Am is the thermal energy change due to gravity and 
hydrodynamics. At is the system time step and A is the cooling function. This equation is 
integrated from to At using many time steps depending on the cooling time scale. Our 
numerical experiments show that sometimes about 10^ time steps is required within one 
dynamic time step At. 



3.3. Mesh Refinement Algorithms 



For an unstructured mesh, cells can be refined arbitrarily. [Vijayan fc Kallinderis (1994) 
discuss strategies for cell division. One can put a node in the center of an edge, in the 
middle of a face, in the middle of a cell, or a combination of all of them. 

In cosmological simulations, we want to resolve forming structures. For N-body 
problems, we use a mass criteria rric to determine mesh refinement. After we interpolate 
particle data to mesh nodes, each node carries a mass. For each face of a cell in the 
unstructured mesh, we put a refining node at the center of the face if the linearly 
interpolated mass at the center is above rric. 

For hydrodynamic problems, following the criteria for galaxy formation in (|Cen fc| 
Ostriker 1993|) , we put a refining node in the middle of a cell if the gas in this cell, (a) 
is contracting, V ■ v < 0, and (b) has a mass greater than the Jean's mass, ms > mj. 
When cooling processes are included we also require the cooling time to be shorter than the 
dynamical time, tcooi < ^dyn- 



4. Code Tests and Performance 

4.1. Testing ID gas-kinetic scheme 

In Figure ^ we show the results of a Lax shock tube test. The initial conditions for 
this test isU = (0.445, 0.311, 8.928) for x < and U = (0.5, 0, 1.4275) for x > 0. The result 
is at t = 0.15. The contact discontinuity is resolved with about 2 cells, the rarefaction shock 
was sharply captured with about two cells and no post-shock oscillation is observed. 



One test that is closely related to structure formation in cosmology is described in Ryu 



et al. (1993)1 with the following initial conditions: p = 0, v{x) = sin(27r2;)/27r, p = 10 and 



periodic boundary condition for < a; < 1. Our results at t = 3 are presented in Figure ^. 
We notice that the BGK gas-kinetic scheme can successfully resolve features within two cells 
without any artificial viscosity or adjustment for the temperature term due to high Mach 
number. Our scheme successfully reproduces the density caustic, the saw-shape velocity 
field, and segmented pressure field with small oscillation. The Mach number involved in this 
test is much higher (~ 10^) than that in usual shock tube tests. This result demonstrates 
that the BGK gas-kinetic scheme is very robust in high Mach number situations. 
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4.2. Cosmological N-body Simulations 

We construct the initial conditions for the Cold Dark Matter model using the 
Zel'dovich approximation (cf. [h^fstathiou et al. 1985| ). Initially, particles are almost 
uniformly distributed. As the system evolves, structure forms due to gravitational 
clustering. More and more massive objects form as time passes by. We show our results 
with a 32^ uniform mesh in Figure 0, and results with mesh refinement from a 32^ uniform 



mesh in Figure |Ty. The final mesh nodes are shown in Figure |TI], which indicates that 
our mesh refinement traces the particle distribution very well. Here, the mesh refinement 
is performed on the faces of tetrahedra. A new node is put in the middle of a face if the 
mass on all the three nodes are above a certain value. In this test case, the critical mass 
is taken to be brrii, where rrii is the mass for each particle. Visually, we can already see 
the great improvement of the resolution with mesh refinement. The two-body correlation 
function ^(r) (see [Peebles 198U| for a definition) for two simulations with and without mesh 



refinement is shown in Figure 0. The resolution improvement of the simulation with mesh 
refinement over the simulation with mesh refinement is well above a factor of 10, while the 
running time between the two simulations for each time step is less than a factor of two. 



4.3. Cosmological N-body + Gas Simulations 



Recently [Frenk et al. (1996]] have proposed a comparison between different cosmological 



hydrodynamic codes. They set up a constrained initial conditions for a CDM model. The 
initial conditions for the results described below are generated from their density field 
using the Zel'dovich approximation ( pfstathiou et al. 1985| ). Readers can compare some 



of our results to others presented in their paper ([Frenk et al. 199"6D with the same initial 



conditions. All the tests shown below are obtained with 32'^ particles and 32^ raw mesh 
with mesh refinement. 



In Figure |13|, we show the density and temperature contours of a slice in the simulation. 
The density field shows some traces of filamentary structure, and the temperature field is 
almost isothermal in high density regions. 



In Figure n , we show the fraction of mass contours of baryonic matter with various 
density and temperature. This figure summarizes the thermal state of the intergalactic 
medium. From this figure we know that most of the baryonic matter stays at the 
average density and in a temperature range of 10'* — 10^ K. Only a small fraction of 
the baryonic matter is in high density regions, and that at high density regions remains 
at high temperatures (~ 10^ K). The material at underdense regions is cold (T <^ 10^ 
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K). The heating and coohng processes change the thermal state of intergalactic medium 
dramatically. One important feature is that the gas at high density region can actually stay 
very cold. The line indicating the heating and cooling balance in Figure |^ indicates the 
thermal states of these high density, low temperature gas. Our simulations with heating 
and cooling processes actually have some fraction of gas in these states. 

5. Summary and Discussion 

In this paper, we described a new cosmological N-body + gas dynamics code based on 
algorithms for an adaptive, unstructured mesh. The novel elements of this code are: (1) the 
mesh construction; (2) solving N-body systems; (3) solving hydrodynamic equations; (4) 
time step estimation and time integration; (5) mesh refinement; and (6) relevant heating 
and cooling processes for primordial gas. 

The mesh construction with periodic boundary conditions is performed using a 
combined Bowyer- Watson algorithm and local transformation algorithm. The initial mesh 
for cosmological simulations is a uniform, staggered mesh. When some refining grids 
are required, new grid points are added to the mesh structure through the incremental 
Bowyer- Watson algorithm, which modifies the previous mesh structure slightly. The 
incorporation of mesh refinement in unstructured mesh gives one way to achieve high spatial 
resolution at relatively low cost. For best results, a good refinement criterion is essential. 
In general, refinement criterion can be derived both on physical and numerical bases. 

Poisson's equation is discretized on an unstructured mesh and solved using conjugate 
gradient method. Particles are interpolated to the mesh nodes using linear interpolation. 
The resulting N-body algorithm is similar to the Particle-Mesh method with CIC 
interpolation. Because each node in an unstructured mesh has more associated cells, the 
new N-body algorithm has slightly higher force resolution than the PM algorithm. 

We solve Euler's equations using finite-volume method. Flux functions are calculated 
using BGK gas-kinetic scheme. The gas-kinetic scheme constructs a time dependent 
distribution in the middle of an edge and calculate flux functions by moments of the 
distribution function. This scheme provides high resolution for shock capturing and is very 
stable for high Mach number flows. 

The new cosmological code solves dark matter and gas dynamics with the same 
resolution. Our tests demonstrate that this code can provide high spatial resolution by 
mesh refinement. We include relevant cooling and heating processes for the primordial gas 
to simulate the evolution of intergalactic medium accurately. 
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Since the internal data structures for an unstructured mesh are homogeneous, 
unstructured mesh codes can be easily parallelized. The difficult part of parallelization is 
the Delaunay triangulation procedure. A parallel version of the code is being developed 
using the portable Message Passing Interface (MPI) library functions. 

The author would like to thank Michael Norman for suggesting the investigation of 
unstructured mesh schemes, and Timothy J. Barth for kindly providing many references 
and some discussions. I greatly appreciate numerous detailed discussions with Kun Xu 
about the gas-kinetic scheme and sharing opinions about computational fluid dynamics. It 
is a great pleasure to thank Lars Hcrnquist and Jeremiah P. Ostriker for their support and 
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For unstructured meshes, we need to determine the relation between a point and a 
simplex (triangle or tetrahedron), and also some geometric quantities, like the volume 
of a simplex and the surface area of its faces. In this appendix, we give the formulas to 
calculated these quantities generally in n-dimensional space. 

Let pi,p2, ■ ■ ■ ,Pn+i be n + 1 distinct points in n-dimensional space. The n-dimensional 
volume of the simplex T with vertices Pi, • ■ ■ ,Pn+i is given by 



Let q be an arbitrary point in n-space. If the simplex T with vertices Pi, • • • ,Pn+i is 
nondegenerate, i.e., if Vol(pi, • • • ,Pn+i) 7^ 0, the numbers, bi, b2, - ■ ■ , bn+i, satisfying 



V bn+l J 

are called bary centric coordinates of a point q relative to simplex T. It can be shown that, 

b, = Det(p„---,g,---p.^O 
Det(pi, • • • ,pk, • ■ -Pn+i) 

Obviously, bk is a linear function of q. Thus, bk indicates the position of q relative to the 
hyperplane containing the facet of simplex T opposite to vertex pk. bk — when q is in 



A. 



Geometric Relations for Unstructured Meshes 




(Al) 




(A2) 
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Hk, bk > when q is on the same side of Hj. from p^, and 6^ < when q is on the opposite 
side of Hk from p^,. As a consequence, we know that the point q is inside simplex T if and 
only if all bj. > 0. The surface area vector Sk of hyperplane Hj. with its direction pointing 
away from pk is 

Sk = -Vo\{T)Vbk{x). (A4) 



It can be shown that the integration of function f{xi, ■ ■ 
simplex T can be expressed as, 

rl rl-bi /-l-bi b„-i 

fix, y, z)<rx = Vo\{T)n\ / dbi I db-,. 



Xn) over the volume of a 
dbnf{xi,---,Xn). (A5) 



If a simplex T is non-degenerate, it has a unique circum-sphere S. Given an arbitrary 
point q in n-space, we can determine if q is inside, outside or on the sphere 5* by the 
following function, 

/ 1 1 ■■• 1 \ 



Det 



InSphere(g, T) 



1 

Wi 

Pi 



Det 




1 

Wn+l 
Pn+1 j 




(A6) 



where w„ = YJi= 



-^p,*- 



InSphere(g, T) > when q is inside InSphere(g, T) = when 



g is on 5 and InSphere(g, T) < when q is outside S. When calculating the value of 
InSphere(g, T), we should be aware of round-off errors ( Parth 1995|) , because the result of 
triangulation could be wrong due to floating point inaccuracies. To avoid the problems 
caused by floating point round-off errors, we calculate the above InSphere(T, q) function 
using the following formula instead. 



Det 



InSphere(T, q) 



P'l P2 



w. 



n+1 
Pn+1 



Det 




(AT) 



Pn+1 



mm 

fc=l,2,---,n+l ' 



where p'^ = Pk — Q and w'f^ = I]"=o^pV "^^^ value of the above function is compared with 
a small number e, instead of 0, to determine the result of the sphere test. For single 
precision floating operations, we use e = 10~^. Our numerical experiments show that the 
above estimate is sufficient to give the correct Delaunay triangulation. Another way to 



avoid the round-off error is use exact redundant expression calculation (see Parth 1995| for 
more discussions). But there are a lot of extra calculations related to the exact redundant 
expression [Fortune fc van Wyk 1993 . 
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B. CoefRcient calculation in BGK formalism 

In this appendix, we will give detailed formulas to calculate the coefficients Aj^, A^, 
A'j^^, A^^, Bp, and other quantities for the BGK scheme described earlier in the paper. 

At the beginning of each time step, we know the fluid state at the two ends of each 
edge and U^. The interpolated fluid state in the middle of the edge, and f/^, which 
are interpolated from left side and right side respectively, can be constructed from the SLIP 
(Symmetric Limited Positive) formulation ( [Jameson 1995| ), which is derived from the local 
extremum diminishing (LED) principle. The constructed fluid state can be expressed as, 

= Ut + \et and = + ^ef , (Bl) 

where = L{AUaj+i/2, ^Ua,j-i/2) is the limited average, L{u,v) is a limiter, and 
AUa,j+i/2 = Ua,j+i — Ua,j- An example is the van Leer limiter, L{u,v) = when 
u ■ V > 0, and L{u,v) = otherwise. The equilibrium distribution functions and are 
constructed from Ul and Lf^ respectively. 

The macroscopic quantities are moments of distributions. We have 
(^) = / + Apippx)dud^. For compatibility, we require U^{x = — |) = U^. After 

some algebra, we get the solutions of the the coefficients A^ and A^, 

< >^ = and < ^^^p A^ = e^. (B2) 

The notation < ■ ■ ■ > will be deflned later in Appendix C. 

The constructed fluid state at the middle of an edge U'-^ is deflned to be, 

U^=[ ^afodud^ = I i^agodud^ + I iJagodudC (B3) 

The equilibrium distribution g^ is constructed from f/*^. Taking the limit t ^ 0, we 
have f/^(x) = / tpad'^i^ + A'^tppx)dud^. For x = ±|, we require f/^(— |) = and 
f/^(+|) = U^. This gives the solutions of coefficients A^^ and A^^. 

The distribution functions f{t,x,u) and g{t,x,u) must be compatible with each other. 
Conservation laws give the following compatibility condition, 

J Mf-9)dud^ = 0. (B4) 

Applying the integrated solution of f(t,x,u) (equation |T^) to the above compatibility 
equation and integrate over the whole time step T, we have, 

1 A^, 



= 1 iJadEdt- ^ if-g) 
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+ 
+ 



2r(l - e-^/-) - T(l + 6"^/^] (< n,V'«V',3 >^>o 4"^+ < w.V^aV'/? >^:<o ^?'') 
-T + r(l - e-^/^)] < ^|J^^Pp >^ 5^ 
1 - e"^/"] (< ^„ >^>o + <^a >u<o) 
-'r-{T + r)e-^/^] < u^^^f, >^>o < ^i^a^P >u<o) (B5) 

These equations give the solutions to the coefficients B/^. 

For convenience, we give the formula to calculate the time integrated flux functions, 



F^{0,t)dt 



dt / uipafit, 0, u)dud$, 



T-t(1 

+T 



) < Ulpa > 



G 



(B6) 



2r(l - e-^/-) - T(l + e"*/-) < m^^^z^^ >^^o < u'^^^p >^^, 



+ 



TV2-Tr + r2(l -e- 



-T/T^ 




5^ < MV'aV'/3 

e-^/") [< >^>o + < uiJa >^<o 
(T + r)e-^/^) [4 < n^V'^.V'^ >^>o < "'V'aV'/? >^<o 



The collision time r can be derived from classical statistical mechanics to be the mean 
free path divided by the rms velocity of atoms. We use the following formula to estimate r 



p 



(B7) 



where AT is the time step and Ci, C2 are constants. We take Ci = 0.01 and C2 = 1 in our 
calculations. The results are not sensitive to the choice of the actual values of Ci, C2. 



C. Velocity moments 

We define the moment of a quantity w of the equilibrium state go as the following. 



where V is the macroscopic velocity, A = p/2p is the gas temperature, and 

D = p(A/7r)(^+^^/^ is the normalization factor. Here, K is the degree of the internal 



variable ^ and is the space dimension Ku 1993. For a polytropic gas, classical statistical 



mechanics gives 7 = {n + 2)/n, where n is the total number of effective degrees of freedom 
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of the molecule: thus a monoatomic gas has n = 3, 7 = 5/3, and a diatomic gas with two 
rotational degrees of freedom has n = 5, 7 = 7/5. For a flow in A^-dimensional space, we 
have K = n- N = ^ - N. 

7-1 



Following the above definition, we obtain the iterative relation. 



< >=V < > 4 



n+l 



n + l 



2A 



< u" >, 



and the following specific values of moments. 



< u > 

<e> 

< > 



1 

V 
K 
2A 

K{K + 2) 

K{K + 2){K + A) 
8A3 



(C2) 

(C3) 
(C4) 

(C5) 
(C6) 
(C7) 



For velocity moments involving integration over half of the velocity space, the above 
iterative relation (equation |(J2D still holds true, except for the following first few moments. 



-erfc(-v^K) 



-WJ: 



^ erfc(-v^K) + , , , 



ierfc(v^V;) 



1 r- e-^^- 
K.-erfc(v^\4) - —= 
^ 2v ttA 



(C8) 
(C9) 
(CIO) 
(Cll) 



The moments of u""il)a and u^ipa'^p can be derived from the moments of m" and We 
explicitly write them out for reference. 



pJ 

( < 7/" > \ 



(C12) 



y i(< > + < X ^2 y 



(C13) 
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For moments of u^ipa and u'^4'a'4^(3 integrating over half velocity space, the above 
expressions are still good except that one must be aware that < m° > may not equal to 1 in 
the above formulas. 



D. Cooling Rates and Reaction Coefficients 

Various cooling processes relevant to primordial baryonic matter have been included in 
the code. Primordial gas is assumed consist essentially entirely of hydrogen and helium. The 
relevant cooling and heating processes are: bremsstrahlung emission, collisional excitation 
of Hq and He^, collisional ionization of Hq, He^ and He~^, radiative recombination of H^, 
He^ and He^~^, and dielectronic recombination of He~^. Ionization equilibrium is assumed 
to determine the fraction in each species. The reaction rates ki (in unit of sec~^) and the 
related cooling rates Aj (in unit of erg cm^ sec~^) are listed below (c.f. Black 19"8T] , |Cen| 



T992i IKatz, Weinberg, fc Hernquist 199q and |Abel et al. 19961) , (Note that T„ = T/10"K 



1. H^ + e^ H+ + 2e 

k, = 5.85 X io-iiri/V^™-^/^(l + T5'/')-i 
Al = 1.27 X io-2iyi/2g-i57809.i/T^^ ^ Tl'^^n^uno 

2. H+ + H^ + -i 

k2 = 8.40 X 10-''T'/^T^'\1 + T°-^)-i 
A2 = 8.70 X 10-''T'/%~'-\l + T^-')~'nenH+ 

3. He^ + e^He+ + e 

^3 = 2.38 X io-iiri/2e-".4/T^^ ^yV2)-i 

A3 = 9.38 X 10-22TV2e-285335.4/T^^ ^ T^/y^rieUHeO 

4. He+ + e /7e° + 7 

^4 = 1.5 X 10-10^-6.353 

A4 = 1.55 X 10-^^T'^-^^^'^n,nHe+ 
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5. //e+ + e ^ //e++ + 2e 

h = 5.68 X 10-12TV2e-631515.0/T^^ ^ ^V2)-l 

As = 4.95 X io-22ri/2e-63i5i5.o/T(i ^ T5'/')-inen^e+ 

6. i/e++ + e ^ He+ + 7 

A;6 = 3.36 X 10-'''T-'/^T^''''{1 + T^-^)'' 
Ae = 3.48 X IQ-^'^T-'/^T^^-^l + T^-')-^n,nHe++ 

7. He+ + 2e^He^ + e 

k, = 1.9 X io-3r-i-5e-™o-o/^(l + 0.3e-™-°/^) 
A7 = 1.24 X io-i3r-i.5g-47oooo.o/T(^ ^ 0.3e-9^°'^° °/^)nen^^e+ 

8. coUisional excitation of 

As = 7.50 X io-i9e-^^«^^«-°/^(l + T^^')nenHO 

9. coUisional excitation of He'^ 

Ag = 5.54 X 10-17T-3-97e-473638.0/T^^ ^ T5^/2)-lnen^e+ 

10. bremsstrahlung 

Aff = 1.42 X 10-%/r^/2ne(nj^+ +nife+ +4nj^e++) 
= l.l + 0.34exp[-(5.5 -logT)V3.0] 

When an ultraviolet (UV) background radiation field is present, photoionization of H^, 
He^ and ife+ is also included. The photoionization rates are defined by 

r^. = r (Di) 

Ji/i riiy 

where J(z/) is the intensity of the UV background, Ui is the threshold frequency and 
crj(z/) is the photoionization cross section for species i. The heating rate associated with 
photoionization is 

H — nnoeno + nHe°^He° + 'nHe+^He+ ergcm~^sec~^, (D2) 
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where 

= r ^^^a,{u){hu - hu,)dv. (D3) 

Besides these radiative coohng processes, we include inverse Compton coohng off the 
microwave background. The inverse Compton coohng rate is given by ( |Ikeuchi fc OstriEei 



19861) , 

Ac = 5.41 X 10"^^neT(l + 2)^ergcm-3sec"^ (D4) 
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Fig. 1. — Relation between Dirichlet tessellation (dashed lines) and Delaunay triangulation 
(solid lines) in 2-dimensional space. The crosses indicate the center of circumcenters of the 
corresponding triangles. 
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Fig. 2. — Example of 2-dimensional unstructured mesh constructed by Delaunay 
triangulation. 1000 nodes are randomly scattered in a periodic box whose borders are 
indicated by dashed lines. 
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Fig. 3. — Illustration of control volume of node with 2-D case on the left and 3-D case on 
the right. On the left panel, the dashed hues, which connects the centers of edges and the 
middle of triangles, indicates the control volume for node 0. On the right panel, only one 
tetrahedron associated with node is shown. The shaded region represents one peice of the 
control volume boundaries. 



-26- 



o 

O 



N=323 



n I I r 



T I I r 



n I I r 



T I I r 







J I I L 



J I I L 



J I I L 



10 20 

Number of Iterations 



J I I L 



30 



Fig. 4. — The convergence rate of our conjugate gradient algorithm. The error decreases 
almost exponentially as the number of iterasions. 
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Fig. 5. — The force between two particles on an unstructured mesh using 32^ nodes. The 
32^ nodes distributed uniformly in a periodic box. The solid line indicates the law. 
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Fig. 6. — Cooling time scale tcooi = E/\A\ with radiative cooling, photoionization and 
inverse Compton cooling for gas at different density and temperature. The UV radiation is 
assumed to be J{v) = 10^^^(i/L/i/) erg/cm^/sec/Hz, and the compton cooling is calculated 
at 2; = 2. 
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Fig. 7. — Lax shock tube test using the 1- dimensional BGK scheme. The result is at time 
t = 0.15. The BGK scheme captures shocks in two cells. 
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Fig. 8. — One dimensional caustic test with BGK gas-kinetic scheme using 100 zones. The 
results are at t = 3. 
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Fig. 9. — Result of a pure N-body unstructured mesh simulation without mesh refinement. 
All the 32^ particles are projected in the X-Y plane. The raw mesh is a 32^ uniform mesh. 
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Fig. 10. — Result of a pure N-body unstructured mesh simulation with mesh refinement. 
All the 32"^ particles are projected in the X-Y plane. The raw mesh is a 32^ uniform mesh. 
Refinement is done by a mass criterion rric — 5 *mi, where is the mass of each particle. 
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Fig. 11. — Distribution of the mesh nodes projected to the X-Y plane. The particle 
distribution of this simulation is shown in Figure 
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Fig. 12. — The two-point correlation function ^(r) for two dark matter only simulations. 
The solid line is for the simulation with mesh refinement, and the dash line is without mesh 
refinement. In both simulations, 32'^ particles are used on a raw mesh of 32^ nodes. The cell 
size for the raw mesh is 1 Mpc/h. 




Fig. 13. — Panel (a) shows the density contours of a shoe in the simulation. The contour 
levels are p/ < p >= 10*^/^ with A: = 0, 1, . . .. Panel (b) shows the temperature contours of 
the same slice with contour levels T = lO^'^^^^'^K with k = 0, 1,2, . . .. 
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Fig. 14. — Mass contour of gas at different density and temperature. The total mass in the 
box is 1. The contour levels are in 10*^^, i = —4, —8. assumed to be 1 



